*Decomposition*
clear

use "Results\tab_B3\Scenario_5\dataset_with_singles_3_and_5_types.dta"


gen ambition=ambition_type_k_5_s

//generate partners' ambition
by couple_id aar, sort: gen am_male = ambition if koen=="1" & relationship==1
by couple_id aar: gen am_female = ambition if koen=="2" & relationship==1
by couple_id aar: egen maxeduc = max(am_male)
by couple_id aar: replace am_male = maxeduc
drop maxeduc
by couple_id aar: egen maxeduc = max(am_female)
by couple_id aar: replace am_female = maxeduc
drop maxeduc


cap drop ambition_type_k_4_s

gen ambition_type_k_4_s=ambition

*Keep only necessary variables
drop aegte_nr faelle_nr civst educ_eika final_educ individual grad_region wage_start_mean_ambition* wage_growth_ambition* wage_start_mean_ambition_s* wage_growth_ambition_s* married cohab 

*Preparation

{
*get 1980 preferences*

*couple matrix

tab am_male am_female if aar==1980 & koen=="1", matcell(am_base) /*It does not matter if we condition on the husband or wife*/

*single men
tab ambition_type_k_4_s if koen=="1" & relationship==0 & aar==1980, matcell(am_single_m_base)

*single women
tab ambition_type_k_4_s if koen=="2" & relationship==0 & aar==1980, matcell(temp)
mat am_single_f_base=temp'

*preference matrix

mat am_pi_base=J(5,5,.)
forvalues i=1(1)5{
	forvalues j=1(1)5{
		mat am_pi_base[`i',`j']=am_base[`i',`j']/((am_single_m_base[`i',1]*am_single_f_base[1,`j'])^(1/2))
	}
}

}

*get 1980 marginals*
mat am_mar_m_base=J(5,1,.)
forvalues i=1(1)5{
mat am_mar_m_base[`i',1]=am_single_m_base[`i',1]+am_base[`i',1]+am_base[`i',2]+am_base[`i',3]+am_base[`i',4]+am_base[`i',5]
}


mat am_mar_f_base=J(1,5,.)
forvalues i=1(1)5{
mat am_mar_f_base[1,`i']=am_single_f_base[1,`i']+am_base[1,`i']+am_base[2,`i']+am_base[3,`i']+am_base[4,`i']+am_base[5,`i']
}




*get later preferences*

forvalues y=2018(1)2018{
tab am_male am_female if aar==`y' & koen=="1", matcell(am_`y') /*It does not matter if we condition on the husband or wife*/


*single men
tab ambition_type_k_4_s if koen=="1" & relationship==0 & aar==`y', matcell(am_single_m_`y')


*single women
tab ambition_type_k_4_s if koen=="2" & relationship==0 & aar==`y', matcell(temp)
mat am_single_f_`y'=temp'


*preference matrix

mat am_pi_`y'=J(5,5,.)
forvalues i=1(1)5{
	forvalues j=1(1)5{
		mat am_pi_`y'[`i',`j']=am_`y'[`i',`j']/((am_single_m_`y'[`i',1]*am_single_f_`y'[1,`j'])^(1/2))
	}
}



*get later marginals*
mat am_mar_m_`y'=J(5,1,.)
forvalues i=1(1)5{
mat am_mar_m_`y'[`i',1]=am_single_m_`y'[`i',1]+am_`y'[`i',1]+am_`y'[`i',2]+am_`y'[`i',3]+am_`y'[`i',4]+am_`y'[`i',5]
}


mat am_mar_f_`y'=J(1,5,.)
forvalues i=1(1)5{
mat am_mar_f_`y'[1,`i']=am_single_f_`y'[1,`i']+am_`y'[1,`i']+am_`y'[2,`i']+am_`y'[3,`i']+am_`y'[4,`i']+am_`y'[5,`i']
}



}

**Inequality**

keep if aar==1980 | aar==2018

gen couple_id_with_s=couple_id
replace couple_id_with_s=pnr if relationship==0
sort couple_id_with_s aar koen

gen individual_earnings=erhvervsindk_13
replace individual_earnings=0 if individual_earnings<0

by couple_id_with_s aar: egen household_earnings=sum(individual_earnings)

**ADDITION**
replace household_earnings=household_earnings/1.5 if relationship==1
**

keep if (koen=="1" & relationship==1) | relationship==0 

*Keep only necessary variables
drop erhvervsindk_13 couple_id couple_id_with_s individual_earnings

**Data**
{

*Data with singles*
foreach y in 1980 2018{
qui ineqdeco household_earnings if aar==`y'

sca gini_`y'=r(gini)
sca p90p50_`y'=r(p90p50)
sca p50p10_`y'=r(p50)/r(p10)
}

frame copy default results
frame change results
gen temp=.
collapse (first) temp, by(aar)
drop temp

gen gini=.
foreach i in 1980 2018{
replace gini=gini_`i' if aar==`i' 
}

gen p90p50=.
foreach i in 1980 2018{
replace p90p50=p90p50_`i' if aar==`i' 
}

gen p50p10=.
foreach i in 1980 2018{
replace p50p10=p50p10_`i' if aar==`i' 
}

frame change default



}



*Fixed household composition*

{
*Solve for 2018 single numbers with 1980 preferences and 1980 marginals (both)*

*ambition types

program nlfaq

syntax varlist(min=1 max=1) [if], at(name)

tempname m1 m2 m3 m4 m5 f1 f2 f3 f4 f5

scalar `m1' = `at'[1,1]
scalar `m2' = `at'[1,2]
scalar `m3' = `at'[1,3]
scalar `m4' = `at'[1,4]
scalar `m5' = `at'[1,5]

scalar `f1' = `at'[1,6]
scalar `f2' = `at'[1,7]
scalar `f3' = `at'[1,8]
scalar `f4' = `at'[1,9]
scalar `f5' = `at'[1,10]

tempvar yh
gen double `yh'= am_pi_base[1,1]*((`m1'*`f1')^(1/2))+am_pi_base[1,2]*((`m1'*`f2')^(1/2))+am_pi_base[1,3]*((`m1'*`f3')^(1/2))+am_pi_base[1,4]*((`m1'*`f4')^(1/2))+am_pi_base[1,5]*((`m1'*`f5')^(1/2))+`m1'-am_mar_m_base[1,1]+1 in 1
replace `yh'=am_pi_base[2,1]*((`m2'*`f1')^(1/2))+am_pi_base[2,2]*((`m2'*`f2')^(1/2))+am_pi_base[2,3]*((`m2'*`f3')^(1/2))+am_pi_base[2,4]*((`m2'*`f4')^(1/2))+am_pi_base[2,5]*((`m2'*`f5')^(1/2))+`m2'-am_mar_m_base[2,1]  in 2
replace `yh'=am_pi_base[3,1]*((`m3'*`f1')^(1/2))+am_pi_base[3,2]*((`m3'*`f2')^(1/2))+am_pi_base[3,3]*((`m3'*`f3')^(1/2))+am_pi_base[3,4]*((`m3'*`f4')^(1/2))+am_pi_base[3,5]*((`m3'*`f5')^(1/2))+`m3'-am_mar_m_base[3,1] in 3
replace `yh'=am_pi_base[4,1]*((`m4'*`f1')^(1/2))+am_pi_base[4,2]*((`m4'*`f2')^(1/2))+am_pi_base[4,3]*((`m4'*`f3')^(1/2))+am_pi_base[4,4]*((`m4'*`f4')^(1/2))+am_pi_base[4,5]*((`m4'*`f5')^(1/2))+`m4'-am_mar_m_base[4,1] in 4
replace `yh'=am_pi_base[5,1]*((`m5'*`f1')^(1/2))+am_pi_base[5,2]*((`m5'*`f2')^(1/2))+am_pi_base[5,3]*((`m5'*`f3')^(1/2))+am_pi_base[5,4]*((`m5'*`f4')^(1/2))+am_pi_base[5,5]*((`m5'*`f5')^(1/2))+`m5'-am_mar_m_base[5,1] in 5

replace `yh'=am_pi_base[1,1]*((`m1'*`f1')^(1/2))+am_pi_base[2,1]*((`m2'*`f1')^(1/2))+am_pi_base[3,1]*((`m3'*`f1')^(1/2))+am_pi_base[4,1]*((`m4'*`f1')^(1/2))+am_pi_base[5,1]*((`m5'*`f1')^(1/2))+`f1'-am_mar_f_base[1,1] in 6
replace `yh'=am_pi_base[1,2]*((`m1'*`f2')^(1/2))+am_pi_base[2,2]*((`m2'*`f2')^(1/2))+am_pi_base[3,2]*((`m3'*`f2')^(1/2))+am_pi_base[4,2]*((`m4'*`f2')^(1/2))+am_pi_base[5,2]*((`m5'*`f2')^(1/2))+`f2'-am_mar_f_base[1,2] in 7
replace `yh'=am_pi_base[1,3]*((`m1'*`f3')^(1/2))+am_pi_base[2,3]*((`m2'*`f3')^(1/2))+am_pi_base[3,3]*((`m3'*`f3')^(1/2))+am_pi_base[4,3]*((`m4'*`f3')^(1/2))+am_pi_base[5,3]*((`m5'*`f3')^(1/2))+`f3'-am_mar_f_base[1,3] in 8
replace `yh'=am_pi_base[1,4]*((`m1'*`f4')^(1/2))+am_pi_base[2,4]*((`m2'*`f4')^(1/2))+am_pi_base[3,4]*((`m3'*`f4')^(1/2))+am_pi_base[4,4]*((`m4'*`f4')^(1/2))+am_pi_base[5,4]*((`m5'*`f4')^(1/2))+`f4'-am_mar_f_base[1,4] in 9
replace `yh'=am_pi_base[1,5]*((`m1'*`f5')^(1/2))+am_pi_base[2,5]*((`m2'*`f5')^(1/2))+am_pi_base[3,5]*((`m3'*`f5')^(1/2))+am_pi_base[4,5]*((`m4'*`f5')^(1/2))+am_pi_base[5,5]*((`m5'*`f5')^(1/2))+`f5'-am_mar_f_base[1,5] in 10



replace `varlist'=`yh'
end


frame create temp
frame change temp

clear
set obs 10

gen y = 0
replace y=1 in 1

nl faq @ y, parameters(m1 m2 m3 m4 m5 f1 f2 f3 f4 f5) initial(m1 1 m2 1 m3 1 m4 1 m5 1 f1 1 f2 1 f3 1 f4 1 f5 1)

*get counter factual singles
mat am_single_m_cf_sort_2018=J(5,1,.)
mat am_single_m_cf_sort_2018[1,1]=[m1]_b[_cons]
mat am_single_m_cf_sort_2018[2,1]=[m2]_b[_cons]
mat am_single_m_cf_sort_2018[3,1]=[m3]_b[_cons]
mat am_single_m_cf_sort_2018[4,1]=[m4]_b[_cons]
mat am_single_m_cf_sort_2018[5,1]=[m5]_b[_cons]


mat am_single_f_cf_sort_2018=J(1,5,.)
mat am_single_f_cf_sort_2018[1,1]=[f1]_b[_cons]
mat am_single_f_cf_sort_2018[1,2]=[f2]_b[_cons]
mat am_single_f_cf_sort_2018[1,3]=[f3]_b[_cons]
mat am_single_f_cf_sort_2018[1,4]=[f4]_b[_cons]
mat am_single_f_cf_sort_2018[1,5]=[f5]_b[_cons]



frame change default
frame drop temp
prog drop nlfaq



*get counter factual couple numbers*

mat am_cf_sort_2018=J(5,5,.)

forvalues i=1(1)5{
	forvalues j=1(1)5{
		mat am_cf_sort_2018[`i',`j']=am_pi_base[`i',`j']*((am_single_m_cf_sort_2018[`i',1]*am_single_f_cf_sort_2018[1,`j'])^(1/2))
	}
}

}


*obtain weights - ambition

frame change results

gen gini_cf_sort=.
gen p90p50_cf_sort=.
gen p50p10_cf_sort=.

gen gini_m_cf_sort=.
gen p90p50_m_cf_sort=.
gen p50p10_m_cf_sort=.

frame change default

forvalues y=2018(1)2018{
	
	
*numerator 1
mat am_num_1_couples=J(5,5,.) /*The prob. of counterfactual given couple cell x,x*/

forvalues i=1(1)5{
	forvalues j=1(1)5{
		mat am_num_1_couples[`i',`j']=am_cf_sort_2018[`i',`j']/(am_cf_sort_2018[`i',`j']+am_`y'[`i',`j'])
	}
}

mat am_num_1_single_m=J(5,1,.) /*The prob. of counterfactual given single men cell x*/

forvalues i=1(1)5{
	mat am_num_1_single_m[`i',1]=am_single_m_cf_sort_2018[`i',1]/(am_single_m_cf_sort_2018[`i',1]+am_single_m_`y'[`i',1])
}

mat am_num_1_single_f=J(1,5,.) /*The prob. of counterfactual given single women cell x*/

forvalues i=1(1)5{
	mat am_num_1_single_f[1,`i']=am_single_f_cf_sort_2018[1,`i']/(am_single_f_cf_sort_2018[1,`i']+am_single_f_`y'[1,`i'])
}

*denominator 1
mat am_de_1_couples=J(5,5,.) /*The prob. of factual given couple cell x,x*/

forvalues i=1(1)5{
	forvalues j=1(1)5{
		mat am_de_1_couples[`i',`j']=am_`y'[`i',`j']/(am_cf_sort_2018[`i',`j']+am_`y'[`i',`j'])
	}
}

mat am_de_1_single_m=J(5,1,.) /*The prob. of factual given single men cell x*/

forvalues i=1(1)5{
	mat am_de_1_single_m[`i',1]=am_single_m_`y'[`i',1]/(am_single_m_cf_sort_2018[`i',1]+am_single_m_`y'[`i',1])
}

mat am_de_1_single_f=J(1,5,.) /*The prob. of counterfactual given single women cell x*/

forvalues i=1(1)5{
	mat am_de_1_single_f[1,`i']=am_single_f_`y'[1,`i']/(am_single_f_cf_sort_2018[1,`i']+am_single_f_`y'[1,`i'])
}

*numerator 2 /*The prob. of being in factual sample*/

sca number_factual=am_`y'[1,1]+am_`y'[1,2]+am_`y'[1,3]+am_`y'[1,4]+am_`y'[1,5] ///
	+am_`y'[2,1]+am_`y'[2,2]+am_`y'[2,3]+am_`y'[2,4]+am_`y'[2,5] ///
	+am_`y'[3,1]+am_`y'[3,2]+am_`y'[3,3]+am_`y'[3,4]+am_`y'[3,5] ///
	+am_`y'[4,1]+am_`y'[4,2]+am_`y'[4,3]+am_`y'[4,4]+am_`y'[4,5] ///
	+am_`y'[5,1]+am_`y'[5,2]+am_`y'[5,3]+am_`y'[5,4]+am_`y'[5,5] ///
	+am_single_m_`y'[1,1]+am_single_m_`y'[2,1]+am_single_m_`y'[3,1]+am_single_m_`y'[4,1]+am_single_m_`y'[5,1] ///
	+am_single_f_`y'[1,1]+am_single_f_`y'[1,2]+am_single_f_`y'[1,3]+am_single_f_`y'[1,4]+am_single_f_`y'[1,5]
	
sca number_counterfactual=am_cf_sort_2018[1,1]+am_cf_sort_2018[1,2]+am_cf_sort_2018[1,3]+am_cf_sort_2018[1,4]+am_cf_sort_2018[1,5] ///
	+am_cf_sort_2018[2,1]+am_cf_sort_2018[2,2]+am_cf_sort_2018[2,3]+am_cf_sort_2018[2,4]+am_cf_sort_2018[2,5] ///
	+am_cf_sort_2018[3,1]+am_cf_sort_2018[3,2]+am_cf_sort_2018[3,3]+am_cf_sort_2018[3,4]+am_cf_sort_2018[3,5] ///
	+am_cf_sort_2018[4,1]+am_cf_sort_2018[4,2]+am_cf_sort_2018[4,3]+am_cf_sort_2018[4,4]+am_cf_sort_2018[4,5] ///
	+am_cf_sort_2018[5,1]+am_cf_sort_2018[5,2]+am_cf_sort_2018[5,3]+am_cf_sort_2018[5,4]+am_cf_sort_2018[5,5] ///
	+am_single_m_cf_sort_2018[1,1]+am_single_m_cf_sort_2018[2,1]+am_single_m_cf_sort_2018[3,1]+am_single_m_cf_sort_2018[4,1]+am_single_m_cf_sort_2018[5,1] ///
	+am_single_f_cf_sort_2018[1,1]+am_single_f_cf_sort_2018[1,2]+am_single_f_cf_sort_2018[1,3]+am_single_f_cf_sort_2018[1,4]+am_single_f_cf_sort_2018[1,5]
	
sca am_num_2=number_factual/(number_factual+number_counterfactual)

*denominator 2 /*The prob. of being in counterfactual sample*/
sca am_de_2=number_counterfactual/(number_factual+number_counterfactual)


*weights
mat am_w_couples=J(5,5,.)

forvalues i=1(1)5{
	forvalues j=1(1)5{
		mat am_w_couples[`i',`j']=(am_num_1_couples[`i',`j']/am_de_1_couples[`i',`j'])*(am_num_2/am_de_2)
	}
}

mat am_w_single_m=J(5,1,.)

forvalues i=1(1)5{
		mat am_w_single_m[`i',1]=(am_num_1_single_m[`i',1]/am_de_1_single_m[`i',1])*(am_num_2/am_de_2)
}

mat am_w_single_f=J(1,5,.)

forvalues i=1(1)5{
		mat am_w_single_f[1,`i']=(am_num_1_single_f[1,`i']/am_de_1_single_f[1,`i'])*(am_num_2/am_de_2)
}

gen weights=.
forvalues i=1(1)5{
	forvalues j=1(1)5{
replace weights=am_w_couples[`i',`j'] if am_male==`i' & am_female==`j' & relationship==1
	}
}

forvalues i=1(1)5{
	replace weights=am_w_single_m[`i',1] if relationship==0 & koen=="1" & ambition_type_k_4_s==`i'
}

forvalues i=1(1)5{
	replace weights=am_w_single_f[1,`i'] if relationship==0 & koen=="2" & ambition_type_k_4_s==`i'
}

replace weights=1 if weights==. /*Then those with missing type are counted same way in factual and counterfactual scenario*/

*Inequality 

qui ineqdeco household_earnings [weight=weights] if aar==`y'

sca gini_cf_sort_`y'=r(gini)
sca p90p50_cf_sort_`y'=r(p90p50)
sca p50p10_cf_sort_`y'=r(p50)/r(p10)



frame change results


replace gini_cf_sort=gini_cf_sort_`y' if aar==`y' 

replace p90p50_cf_sort=p90p50_cf_sort_`y' if aar==`y' 

replace p50p10_cf_sort=p50p10_cf_sort_`y' if aar==`y' 


frame change default



drop weights
}



*SAVE THE RESULTS*
frame change results

save "Results\tab_B3\Scenario_6\Decomposition\results_eq.dta", replace


